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ABSTRACT 


In all practical application of plasma, the sheath formed close to wall material plays an crucial 
role in determining the property of over all plasma wall transition. Mostly collisionless sheath 
were studied using fluid model and kinetic model but it is the first time we try to study 
collisional sheath using Kinetic Trajectory Simulation model by introducing general collision 
term explicit dependence on a constant a, which determines the degree of collision. In this 
model, the distribution function of the particles species are directly calculated by solving 
the related collisionless and time-independent kinetic equation along the respective particles 
trajectories. We vary the value of alpha from 0 which is for collisionless case to 1 and study 
the different sheath parameters. But increasing a after certain value there is no significant 
change in sheath parameters due to saturation in sheath. The result obtained is expected to 
be useful in controlling flow of plasma particles at the wall. Our study has crucial importance 


in material processing, plasma etching and for confinement of plasma in fusion devices. 


Chapter 1 


Introduction 
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1.1 Plasma 


The term plasma was introduced by Langmuir in 1928 to describe the state of matter in the 
positive column of glow discharge tubes [1]. The word plasma means something molded or 
fabricated ( comes from Greek word HAaopa ). The existence of plasma as the fourth state 
of matter was first identified by Sir W. Crooks in 1879 ( he called it radiant matter ) [2]. 
If temperature of gas increased beyond a certain limit , it does not remain a gas, it enter a 
region where the thermal energy of it’s constituent particle is so great that the electrostatic 
forces which binds electron to atomic nuclei are overcome. Instead of hot gas composed of 
electrically neutral atoms, we have mixed population of charged and neutral particles. This 
is a plasma, and it is neither solid, liquid nor gas. But any ionized gas can not be plasma 
[3], there is always some small degree of ionization in any gas given by Shahas Equation. 
The amount of ionization to be expected in a gas in thermal equilibrium is given by Saha 
equation: 
ae 2.4 x 192? visor (1.1) 
n” n 
where n'’ and n” are respectively, the number density of ionized atoms and neutral atoms, T 
is the gas temperature in Kelvin, kg is Boltzmann’s constant, U; is the ionization energy of 
the gas. For ordinary gas at room temperature, we can take n” = 3 x 107° m-3, T = 300° 
K, and U; = 14.5 eV (for Nitrogen). The fractional ionization is low: 

ie ~~ 10°” (1.2) 

nr 
As the temperature is raised, the degree of ionization remains low until U; is only a few 
times kgT. Then n‘/n” rise abruptly, and the gas is in a plasma state. Further increase in 
temperature makes n” less than n', and the plasma eventually becomes fully ionized. 
Thus plasma is quasineutral gas that shows collective behavior. Quasineutrilaty means elec- 
tron density and ion density are nearly equal ( n; © n. © n, where n is density of plasma 


) but not so neutral that all the interesting electromagnetic vanish. Collective behavior of 
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plasma implies that the dynamics of plasma particles depends not only an the local condition 
but also on the state of plasma far away from the point of interest [3]. 


Any ionized gas can be called plasma if it fulfill the flowing criterion 


1) Ap << LE (Debye shelding). 

2) The number of particle in Debye sphere is very large, Np >> 1. 

3) wpT > 1, where wy is plasma frequency and 7 is mean time between the collision with 
neutral atoms. 

A fundamental property of plasma is its ability to shield out any external electric potential 
applied to it. Before external applied potential to plasma, the plasma potential ¢ = 0 every 
where because ion and electron densities ( -ve and +ve charge density ) were spatially uni- 
form and equal. The particle having same polarity as the test charge ( used to apply external 
potential ) will be repelled where as particle of opposite polarity will be attracted. In the way 
the test charge surrounded by the opposite charge to that of test charge to alleviate it’s effect 
on remaining bulk of plasma. i.e external potential shielded in a small. This phenomenon is 
called Debye shielding. This region is generally considered as spherical in shape being test 
particle at centre, known as Debye sphere and radius of this sphere is called Debye length. 
As shown in figure 1.1. 

The typical expression for the electron Debye length with electron temperature T° and com- 
mon density n° is given by 


" €okpIé 
\S = 4/ ae (1.3) 


The picture of Debye shielding |cf. Fig. 1.1] that we have given is valid only if there are 


enough particles in the charge cloud. Clearly, if there are only one or two particles in the 


sheath region, Debye shielding would not be a statistically valid concept. 


4 


As ion are massive than electron can be considered as uniform back ground in plasma ( 
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® e® @e 
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Figure 1.1: Debye shielding. 


plasma approximation ). If the electron in plasma are displaced uniform back ground of ion, 
a electric field will be built up in such a way as to restore the neutrality of the plasma by 
pulling the electron back to their equilibrium position. Because of their inertia, the electron 
will be overshoot and oscillate around equilibrium position with a characteristic frequency 
known as plasma frequency. The expression for plasma frequency is given by 


Noe? 


(1.5) 


WwW = 
. EgMe 


Where n, is plasma density,e and m, respectively charge and mass of electron. 

The third condition imply motion of plasma particles is controlled by electromagnetic forces 
rather than by ordinary hydrodynamic forces. 

The study of plasma is crucial due to its application in various field like astro physics, 
material processing,surface treatment, fusion energy, semiconductor device etc. One of vital 
source of energy is sun in which nuclear fusion reaction takes place and the energy produces 
continuously same as in stars [5]. Various methods for the terrestrial realization of fusion 
energy have been explored so far but confinement of plasma is one the challenging problem in 
fusion. In such fusion devices where the plasma comes into contact with a wall, the plasma 


region near the wall plays an important role in determining the overall plasma properties 
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which is known as the “sheath”. 


1.2 Sheath 


In the production of plasma energy supplied to gas in order ionize is not problem but problem 
is that how to confine it for a certain time. For any practical use it must has to come in 
contact with the material. When plasma comes in contact with suitable material due to high 
mobility of electron ( due to lighter mass than ion ) they hits the wall charging the material 
surface negative respect to bulk of plasma. As the negative potential of wall increases ions 
are attracted and it repels the part of electron. An equilibrium is finally reached when the 
potential difference is few times the electron temperature forming a positive space charge 
region is in between wall and bulk plasma known as sheath. Simply it can be stated that the 
sheath region is formed due to Debye shielding so scale of sheath will be few Debye length Xp. 
The function of a sheath is to shield quasineutral plasma from the potential distortion caused 
by an absorbing negative wall. Due to the shielding effect of plasma, the negative potential 
at wall have stronf effect only in the sheath region which extended only to few Debye length 
away from the wall. With in this region the plasma is significantly non-neutral how ever 
becoming quasineutral at the sheath edge ( also called “sheath entrance ( SE )” ). As we 
move very close the wall, the potential falls off rapidly. Due to decrease in the potential the 
electric field is strong and the motion of the particles is dominated by electric force rather 
than magnetic force. The sheath structure is responsible for the flow of the particles and the 
energy towards the wall and may also affect the bulk-plasma behavior [6, 7]. The sheath is 
dominated by large gradient of the electric field compared to that of presheath region. The 


Fig. 1.2 shows the schematic potential variation in front of a negative wall. 
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Presheath Bulk 
Sheath Plasma 
——©® 
Fast 


Figure 1.2: Schematic potential variation in front of a negative wall. 


1.3. Presheath 


A quasineutral region in between sheath and bulk plasma having length scale greater than 
sheath and less than bulk plasma formed [cf. Fig. 1.1] due to presence of weak electric field, 
which is due to not perfect shielding of wall potential, responsible for accelerating the ion to 


a velocity that satisfy Bohm criterion. The Bohm criteria [8] is 


1 


DS a (1.6) 
where () denotes averaging over the ion distribution function and 
k se + fA 
en ei pe) (1.7) 
m* 


is the ion-acoustic velocity defined at the presheath side of the sheath edge. Here y’ and ¥° 
are the ion and electron polytropic constants respectively and DS. and Ty, are the ion and 


electron temperatures at the presheath side of the sheath edge respectively [6, 7]. 


CHAPTER 1. INTRODUCTION t 


Thus to form sheath ion must enters the sheath region with minimun or greater than ion 


acoustic velocity ( C, ). 


1.4 Collision in plasma 


Study of collision process between micro particles of a system like fluid and plasma has great 
importance became such systems are usually associated with the collisions. By study the 
collision, we can investigate the nature of the force existing between them. Not only in- 
vestigate the nature of force, we can study different phenomenon in plasma. The collisions 
are also responsible for bringing the system about the ultimate Equilibrium situation where 
the Maxwellian distribution prominent. The collisions between the particles may be elastic 
or inelastic. Elastic collisions, also called soft collisions or grazing collisions are collision in 
which there is no change in the internal energy of the colliding particles, these collisions are 
studied classically. In inelastic collision there is change in the internal energy of the collid- 
ing particles and they are studied quantum mechanically. It is noted that whether collision 
is elastic or inelastic there is always transfer of momentum between them. In the collision 
of neutral particles momentum transfer occurs only impulsively where as in the collision 
of charged particles momentum is transferred continuously through grazing collisions. In 
the classical theory of collision, particles are taken as rigid and spherically symmetric point 
masses, having no spin and so having only three degree of freedom due to only translational 
motion.The average of distances travelled between successive collisions is called the mean 
free path. Mean free path is taken as the average of distances travelled during a sequence 
of small angle collisions. The average time that spends by particle between two successive 
collision is called the collision time. 

In weakly ionized plasma collision infraction between charged particles and neutral particles 
and between neutral particles and short range infraction. But the charged particles them- 


selves infract mutually with a long range coulomb infraction. So, in weakly ionized plasma, 
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we consider both short and long range infraction. In fully ionized plasma, particles mutually 
infract with long range Coulomb force. Due to such long range infraction particles undergoes 


continuous deflection during their motion [9]. 


1.5 Collision in sheath 


The collision in sheath not similar to collision in plasma because sheath is non-neutral due to 
high density of ions than electrons, sheath itself alone can not be called plasma. The collision 
in sheath region result of high density of ions, thus it is close to sufficient to account collision 


between ions. 


1.6 Scope and applications 


The collisonal effect in sheath region have studied with help of Kinetic Trajectory Simulation 
(KTS) [16] with slight changes to include collision; as this model is expected to give more 
accurate result. In all practical application of plasma, it must come in contact with material 
where a sheath layer formed between material wall and bulk plasma. So it is crucial to have 
knowledge of sheath, different sheath parameter and different sheath phenomenon to work 
in field of plasma. 

The plasma physics applied in various field like in nuclear physics, astro physics, surface 
treatment of material, plasma diagnosis etc. 

The dissertation work is organized as introduction is given in chapter 1, in chapter 2, reviews 
of some related work are discussed. In chapter 3, the principle of KTS is discussed. In 
chapter 4, we describe the /d1v plasma sheath model with boundary conditions, discuss the 


kinetic Bohm-Chodura criterion and the coupling of presheath and sheath. In chapter 5, we 
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discuss the numerical methods for the plasma sheath simulation where we have also described 
how our solution is iterated, starting from given boundary conditions to reach the final self- 
consistent state. In chapter 6, we have presented the results obtained with discussion and 


conclusion and numerical program (MATLAB) files are listed in the Appendix. 


Chapter 2 


Review of Plasma Sheath 


CHAPTER 2. REVIEW OF PLASMA SHEATH Et 


Studies on confinement of plasma have always been a priority of interest for all practical 
plasma devices, specially for the magnetic confinement fusion device, where the plasma is 
confined in a vacuum chamber of a finite size. The study of plasma at the wall of the cham- 
ber has not yet been fully understood. However, several theoretical and practical works have 


been done. Some of which are reviewed below. 


T . E. Sheridan and J. Goree investigated collisional sheath model [10], the effects of ion 
collisionality on the plasma sheath are revealed by a two-fluid model. The ion-neutral colli- 
sion cross section is modeled using a power law dependence on ion energy. Exact numerical 
solutions of the model are used to determine the collisional dependence of the sheath width 
and the ion impact energy at the wall. Approximate analytical solutions appropriate for the 
collisionless and collisionally dominated regimes are derived. These approximate solutions 
are used to find the amount of collisionality at the center of the transition regime separating 
the collisionless and collisional regimes. Rx- the constant ion mean-free-path case, the center 
of the transition regime for the sheath width is at a sheath width of five mean-free paths. 
The center of the transition regime for the ion impact energy is at a sheath width of about 


one-half of a mean-free path. 


S.Farhad Masoudi, Shadi S$. Esmaeili, Shima Jazavandi study the ion dynamics in plasma 
sheath under the effect of EF x B and collision force using fluid model.They investigated the 
velocity [11], kinetic energy and density distribution of ion in collisional and collisionless 
magnetic plasma sheath. Considering an external magnetic field, the ion moment under the 
effect of magnetic, electric and collisional force has been analyzed. It is shown that sheath 
thickness decreases by increasinf the magnetic field strength and effect of collision on ion 


characteristics illustrated. 


I. Langmuir was the first study of plasma by the basic features of the plasma sheath transi- 


tion [7]. He studied the interaction of electron and positive ion space charges in sheaths. He 
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found that the potential distribution in the plasma, determines the motion of the ions and 
thus fixes the rate at which the ions arrive at sheath. In his investigations, the essence of 
Bohm criterion was used in an implicit form. In 1949, Bohm formulated and interpreted the 


sheath formation condition explicitly [8, 19]. 


H.-B. Valentinia study that in collisional in slightly ionized plasmas, i.e. for many typi- 
cal low pressure discharges, at the sheath edge the drift velocity of the plasma can be smaller 
than the ambipolar sound speed of the ions. This basic result demonstrated using simple 
analytical expressions. Under these conditions the Bohm criterion, in its usual formulation, is 
a sufficient condition for the existence of a positive boundary sheath, but it is not a necessary 


condition [13]. 


kK. U. Riemann, F. Hermann and H. B. Valentini attempted to generalize the Bohm cri- 
terion for collisional case [?, ?, ?]. These attempts were based on an arbitrary definition of 
presheath and sheath regions. Riemann established a well-known and widely used subdivi- 
sion of the boundary layer of a collision dominated plasma into a collisionless sheath and 
quasineutral presheath. This subdivision is strictly valid only in the case when au — 0, 
where Xp is the electron Debye length and A is the ion mean free path [5]. The sheath edge 
separating the presheath and sheath is defined by the Bohm criterion. Riemann has used a 
simple fluid model of ions to account for collisions and space charges in boundary layer to 


clear the controversial statements on the validity of plasma sheath concept and on the role 


of Bohm criterion for finite value of AD. 


M. Goswami and H. Ramachandran illustrated a 1-D self consistent plasma column that 
includes both the sheath and the presheath [14]. They derived an analytic solution for inter- 
esting limit where the collision operator is independent of flow of velocities. They claimed 
that the most conspicuous feature of the boundary potential profile is the plasma sheath 


that is normally formed at plasma boundary to balance electron and ion flux loss. Electric 
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fields at sheath are normally much stronger than bulk plasma electric fields. In collisionless 
and weakly collisional systems, sheath scale lengths depend on Debye length, Ap, and are 


normally much smaller than the plasma dimension, L. 


Rui Rosa demonstrated a simple model of a collisional plasma sheath, with ionizing col- 
lisions included and in contact with a negative metal wall. A broad range of plasma sheath 
solutions is obtained and their properties are discussed. The characteristic of the abnormal 
glow discharge as well as the characteristic of the electrostatic probe, in the ionsaturation 
region, are incidentally recovered. This then applied to the study of current multiplication 
in electrical discharges. The properties of a plasma current multiplier are predicted, both in 


the low and high pressure regimes [15]. 


R. Khanal contributed in the development and understanding of the sheath by developing 
a kinetic model to study the space charge sheath adjacent to an absorbing wall. He cou- 
pled the kinetic sheath solution to a two-fluid presheath one. He described extensively the 
Kinetic Trajectory Simulation (KTS) model and how to obtain the final self consistent time- 
independent state. He devised and demonstrated this model only for 1d1v, time-independent 
and collisionless cases and applied this model to the special case of a /d1v single emitter 
plasma diode for the purpose of testing and comparison. He used the KTS model to study a 
space charge sheath adjacent to an absorbing wall and bordering on a presheath described in 
a two-fluid model. The solution of related presheath-sheath coupling problem is developed, 
evaluated and discussed in detail. A comparison of the sheath model developed by him and 
the ‘standard’ one (having Boltzmann distribution electrons) shows qualitative agreement 


but also some discrepancies due to the different electron distributions [16]. 


V. Godyak suggested that the most prominent feature of the boundary potential profile 
is the plasma sheath that is normally formed at the plasma boundary to balance electron 


and ion flux [18]. According to him, Sheath electric field are much stronger than the bulk 


CHAPTER 2. REVIEW OF PLASMA SHEATH 14 


plasma electric field and in collisionless and weakly collisional systems, sheath scale lengths 


depends on Debye length and are normally much smaller than the dimensional of plasma. 


R. N Franklin study the conventional equations describing plasma and sheath are exam- 
ined over a full range of collisionality to determine the influence of collisionality on the ion 
equation of motionthat is, when the ion motion is essentially inertial and when it is colli- 
sional. When it is inertial in the sheath the Bohm criterion has been satisfied, while when 
the sheath is collisional the ions remain in equilibrium with the field all the way to the wall 


[17]. 


S. Farhad Masoudi investigated collisional sheath using external magnetic field [12], by con- 
sidering collision between ions and neutral gas atoms under various pressure. He showed that 
ion density increase as collision effect than collision less case. Also in sheath near to wall the 


effect of magnetic field dominated by collision. 


Chapter 3 


Principle of Kinetic Trajectory 
Simulation (KTS) 
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3.1 About KTS model 


In the Kinetic Trajectory Simulation (KTS) the velocity distribution function of particle 
species involved are directly calculated by solving the related kinetic equations along the 
respective collisionless particle trajectory. In order to obtain the distribution function at any 
point of the phase-space we trace the related trajectories of phase-space where the distribution 
function is given. Here we assume the electron and ion velocity distribution function at the 
sheath edge to be cut-off Maxwellian. 

KTS is an iterative method for numerically calculating self consistent, time-independent 
kinetic plasma states in some given bounded spatial region. The plasma states are generally 


characterized by, 


e the velocity distribution function, f(x, v), 


the electric field, E(x), 


the magnetic field, B(x), 


collision term, C(x) 


e given boundary conditions. 


3.2 Basic concepts of kinetic theory 


At any given time, each particle has a specific position and velocity. We can therefore charac- 
terize the instantaneous configuration of a large number of particles by specifying the density 
of particles at each point r,v in phase-space. The function prescribing the instantaneous den- 
sity of particles in phase-space is called the distribution function and is denoted by f(r, v,t). 
Thus, f(r,v,t)drdv is the number of particles at time t having positions in the range be- 
tween r and r+dr and velocities in the range between v and v+ dv. As time progresses, the 


particle motion and acceleration causes the number of particles in these r and v ranges to 


CHAPTER 3. PRINCIPLE OF KINETIC TRAJECTORY SIMULATION (KTS) 17 


change and so f will change. This temporal evolution of f gives a description of the system 
more detailed than a fluid description, but less detailed than following the trajectory of each 
individual particle. Using the evolution of f to characterize the system does not keep track of 
the trajectories of individual particles, but rather characterizes classes of particles having the 


same r,v [6]. In the general case of time-dependent, collisional kinetic theory, the species-s 


Vv 


e a es 
e « rm ee eo? 214 
- im = a = : “ 5 : x 


Figure 3.1: A box with in phase space having width dx and height dv in /d1v. 


velocity distribution function satisfies the kinetic equation 


dfs — (a 0. .» OV 
= (Ftv gta ss en 
ee (3.2) 
with 
a’(r,v,t) = — [E(r,t) + v x B(r, 8). (3.3) 


Here, the macroscopic (i.e. space locally averaged) electric and magnetic fields are E(r, t) 
and B(r,t), the macroscopic acceleration of the species-s particles (i.e. its acceleration in 
these fields) is a*(r,v,t), and C® is the species-‘s’ collision term. The time derivative 


d_(0, a, 8 
ae(gty pat 5) (3.4) 


is the “Lagrangian” or total time derivative along the species ‘s’ trajectory. 
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3.2.1 For collisionless case 


For collisionless cases the kinetic equation Eq. (4.2) takes the well known form of “Vlasov 


equation” 
a) 0 OO Ve ie. 
(Gt gta psa ee) 
Le., 
d os 
a 
fo = J. = constant: (3.6) 


where f;, is startiting distribution for species- s. 

This means that the velocity distribution function is same as starting or constant for an 
observer moving along a collisionless trajectory. Hence, the distribution function at every 
point along the trajectory can be obtained if its value at one point is known. Here we assume 


the boundary distribution function is given. 


3.2.2 For collisional case 


For collisional case the kinetic Eq. (4.2) takes the form 


a a eee 

(Gt Rta zt = is) 
L.€., 

d s — s 

“aS 

oS a (3.8) 


where f%, is startiting distribution for species- s and [°* is collision integral given by 


= / C*dt (3.9) 


This means that the velocity distribution function is not constant but changes for an observer 


moving along a collisionless trajectory. Hence, the distribution function at every point along 
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the trajectory can be obtained only if its value at one point and the effect of collision is 


known. Here we assume the boundary distribution function and collision term C’® is given. 


3.3. Basic equations 


The following basic equations correspond to one-dimensional, time-independent, collisional, 


electrostatic problems: 


(a) For electrons (¢q = —e), C® = 0 the velocity distribution function satisfies the time- 


independent Vlasov equation in differential form as 


dfe | 0 € OF ccd 
a Vea — 2, (E(t) + x Bir): f(r, v) (3.10) 


Ce 


= 0 (3.11) 
and in trajectory integrated forms as: 
f(z, v) = ape (Toe Vee) ’ (3.12) 


where f< is the “starting distribution”, r¢, is the starting configuration and vé§, is the starting 


velocity of electrons. The electron equations of motion are 


GEe eee (3.13) 


qe (3.14) 
And 

dv° Z 

i (3.15) 


with component of accelerations 


dus e 
HO te (3.16) 


CHAPTER 3. PRINCIPLE OF KINETIC TRAJECTORY SIMULATION (KTS) 20 


The macroscopic acceleration is 


ae) = [E(r) + v x B(r)| (3.17) 


e€ 
me 
with its components 


a’ = — [E(x) + (v x B(r)),], (3.18) 


(b) Similarly, for the singly charged ions the velocity distribution functions satisfy the time 
independent Vlasov equations in differential form as 


a = |v +5 (@) +vx BO) S| few) (3.19) 


C' (3.20) 


and in trajectory integrated forms as: 


f'(r,v) = i (rey. Vin) a3 I'(r,v) (3.21) 


a 


‘is the starting configuration, v’, is the starting 


where f%, is the “starting distribution”, r 


velocity and J’ collision integral of ions. The ion equations of motion are 


dr’ 
er 22 
ama (3.22) 


with components 


die - 4 

— = 0h, (3.23) 
And 

dv' F 

oo (3.24) 


tq (3.25) 
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The macroscopic acceleration is 


a'(r) E(r) + v x B(r)] (3.26) 


e€ 
~ ine | 
with its components 


ai. = — [E(x) +(v x B(r)),.J, (3.27) 


(c) The electric field is given by 


do(r) 
E —— 3.28 
(r) = -S (3.28) 
where the electrostatic potential ¢(r) is to be found from Poisson’s equation 
2 
dor) _ p(t) (3.29) 


dr? Eo 


The space charge density is defined as 
p(t) = >> ¢n*(r;) (3.30) 
with the electron and ion densities given as 


WAP) = [ du f*(r,v); S= (6:1): (3.31) 


(oe) 


3.4 Boundary conditions 


In the present case, the boundary potential, ¢(r,:), and boundary injection distribution func- 
tion, f% (1st, Vst), are assumed to be given. We choose f% (rst, Vst) such that the distribution 
function describes the physics we are interested in. Particularly, we assume the electron and 
ion velocity distribution functions to be cut-off Maxwellian in order to fulfill the most im- 
portant requirements for presheath-sheath transition, i.e., quasi-neutrality, the sheath edge 
singularity condition, continuity of the first three moments of each species and the kinetic 
Bohm criterion. This model helps to obtain the distribution function at any point (r,v) of 


the phase space if the distribution function at any point in the simulation region is given. 


Chapter 4 


The Plasma Sheath Model 
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4.1 Description of the model 


Our model consist of non-magnetized electrostatic plasma collisional sheath bounded by 


x =0 and x = Las shown in Fig. 4.1. 


Wall Sheath 
Entrance 


Figure 4.1: Plasma sheath model. 


Plasma consists of only electrons and singly charged ions. The right hand boundary (a = L) 
is specified as the “sheath entrance” which separates the non neutral, collision sheath region 
(x < L) from the quasineutral, presheath region (x > L), whereas the left hand boundary 


(x = 0) represents an absorbing wall. We assume collision term such that collision integral 
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becomes 


C.L 


=a 
Xi 


(4.1) 


where 0 < a < 1 represents the degree of collision. 

The value of a@ very over the range to study the sheath in different degree of collision. The 
value of a represents the degree of collision, i.e. how frequent the collision occur. Lets take 
a = 0.01 which implies that if a particle encounters 100 other particles it undergoes a single 
collision, @ = 1 corresponds to the case where every encounter is collision . The a = 0 
represents the collisionless case and 0 < a < 1 represents the collisional casses. The a = 1 
represents fully collisional sheath which is not practically possible. The sheath is almost 


collisionless so it is suitable to take smaller value of a close to 0. 


4.2 Boundary conditions 


To solve the set of equations compiled in Section 3.3, we need to know the boundary condi- 
tions for the velocity distribution functions (particle boundary condition) and the potential 
(field boundary condition) at the two boundaries of the simulation region. The boundary 
potentials, ¢(z = 0) and ¢(a = L), the boundary injection distribution function, f*(L, v) 
and the distribution function at the wall, f*(0,v), are assumed to be given. Hence, they 
must be specified before the iteration is started and are kept constant throughout the entire 


simulation. 


(a) Particle boundary conditions 


We assume here the plasma particle enters the simulation region from the right hand bound- 
ary with cut off Maxwellian velocity distributions functions, the left hand boundary does 
not emit any particle and that both boundaries are perfectly absorbing. Due to this, the 


distribution function satisfies the following boundary conditions: 
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at the left hand boundary, 
eS 2 Urea e= (e;1): (4.2) 


At the right hand boundary, 


mov" 
(ede = A® ——___ 
f(z ,v <0) erp | | 
2 
= Avexp |— = : (4.3) 
Vir 
and 
ria =L,v <0) — A m'(v-vip)? O(v! => ) 
= L. exp Dhar Vet — Vx 
2p 
= A'exp |— — O(vi, — ve), (4.4) 
Vif 


where T° and T" are the electron and ion formal temperatures respectively. The thermal 


velocity of species‘s’ is given by 


, kp? 
Vif = ae 5 (4.5) 


where v’,, is the ion “Maxwellian-maximum” velocity at x = L and vi, (v’,; < 0) is the 


ion cut-off velocity at x = L. We assume the particles enter the injection plane with half 
Maxwellian and our potential profile is monotonically decreasing towards the left hand bound- 


ary. Together with these two conditions the total electron velocity distribution function is 


given by 
vw? ed(a 
Fae = Lvs 0) = Aver uP ae O (ue, — Uz), (4.6) 
t 
where 
9 _ 
Ver =: 2e|o(x) = bo] (4.7) 


is the electron cut-off velocity at point x. For « = L, we have 


“i? 
UE,” 


f° (a = L,v) = A®exp O(ve, — Uz), (4.8) 
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—2ed, 
jaa aes (4.9) 


me 


Similarly, for the ions from Eq. (4.4) and using the fact that our potential profile is monoton- 
ically decreasing towards the wall (left hand boundary) and ions enter the sheath entrance 
with half Maxwellian velocity distribution function gives the total velocity distribution func- 
tion for ions at x = L is 
oe 

file =L,v) = A’exp | — (| O(vi, — v2). (4.10) 
In the right hand side of the Eq. (4.8) and Eq. (4.10), there are seven parameters A°, T¢, u,, 
A’, T;, Ui, and vt; which must be specified according to the physical situation considered. 


Now the particle density and other physical parameters at « = L are given by 


= i du f*(L, v). (4.11) 


The average or fluid velocity, 


s —s s 
tS aS, ye 


/ “aoup Ea (4.12) 


(oe) 


and effective temperature, 


7 2 mi(v—u)? 
Leh = k < D >L (4.13) 


The electron density at x = L is calculated with Eq. (4.11), (4.8) 


ny io duf*(L,v) 


(oe) 


I| 
aN 
av) 
‘aes 
8 8 
Q 
e 
fav) 
8 
Ss 
| 
aN 
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8 eo 
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using (=) = €, we get 
tf 


2 2 UE 
a Ae qh/2 UE; a 
9 D) 


ee ee ee (4.15) 
We define “erf’ as “Error function”, which is given by 


erie) = = if dé exp (-€7). (4.16) 


Now, from the velocity distribution function Eq. (4.10) and Eq. (4.11) we get ion density as 


ni. iz du f'(L,v) 


(oe) 


I| 
C—— 
+ 

8 

Q, 

e 

a 

o 

8 

3 
| 
, 
< 
< | | 
a] < 
3 
w& 
Ssh 
bo 

oO) 

— 

ie) 

w& 
| 

e 

a 


I| 
a 
es 
8g 
Q. 
eS 
| 
i 
S 
8 
SS | 
Sle 
3° 
sp 
NL 
bo 


wh Se ae (aot — al (4.17) 
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where C* = 1+erf(ri,) and 7, = Gane . We can also define the “complementary error 
function” as 


erfeday=1l—erf (x). (4.18) 


The average electron velocity is calculated with help of Eq. (4.8), (4.14), (4.12) 


L, eee 
= — | duv f°(L,v) 
DL J-oo 


_ uf, DE 
where, 
ee Ver 
D(T;, 60) = exp|-—>| 
ego 
= erp[-—2] (4.20) 
kT 
The ion average velocity is given by Eq. (4.10), (4.17), (4.12) 
1 ee eas 
uy = => duv f'(L, v) 
ny —oo 
i YD" 
cod UmL = m/2C% (4.21) 


DU(T}) = exp[-(“ ney) 


= exp[—t?| (4.22) 
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The electron effective temperature is given by Eq. (4.8), (4.14), (4.13), (4.19) 


e mé a e\2 re 
Lt. = me du(v — uf)" f°(L, v) 
2u’, D° By ks 
= Te\1l—- c — 
fl ie) Ox a De 
where, 
meve 2 
( 
2k 


The ion effective temperature is given by Eq. (4.10), (4.17), (4.13), (4.21) 


i ae aie i \2 pi 
Tepe Re es du(v — uh)*f*(L, v) 
: ey Oem Ae 
ae Ei Te. ae ~ (aa) | 
where, 
mivi? 
Ti & tf 
f 2k 


(b) Field boundary conditions 


29 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


The potential profile at « = L is chosen as zero, whereas the one at x = 0 is fixed to a 


negative constant value, i.e. 


o(x = 0) = 5 (constant) < 0 


o(e@=L) = o(L) 
= 0. 


(4.27) 


(4.28) 


We restrict ourselves to the potential distribution which decreases monotonically from 7 = L 


to « = 0 such that the electric field is always negative as shown in Fig. 4.2. 
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Presheath Plasma 
(Collisional) 


Wall 


(%,) 


Figure 4.2: Example of potential profile decreasing monotonically from x = L towards x = 0. 


4.3. Bohm criterion 


In order to form sheath the ions at sheath entrance must have a velocity greater or equal 
ion acoustic velocity which was shown by Bohm and known as “Bohm Criterion”. Here we 
derive Bohm criterion [7, 19], which is satisfied by injected ion velocity distribution function, 
so that formation of sheath is possible, whereas the electron distribution is assumed to satisfy 


full Maxwellian. 


We have the Poisson’s equation 


dx? oa 


do p(@) (4.29) 


Expanding the above equation in Taylor’ series about x ~ L, we get 


do = 1 dp & dp | a 
ales = pte) a SI. oD al. Cie (4.30) 


As we have assumed at the sheath entrance the plasma is quasineutral, i.e. p(x,) = 0, 


nig =n, - Hence, using the quasineutrality condition and neglecting second and higher order 
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terms from Eq. (4.30), we obtain 


dd _ o | dp 5 
ale "es iB] -° (4.31) 


But we have taken potential profile monotonically decreasing towards the wall i.e., of the 


form ¢(a < L) < 0. Hence the Eq. (4.31) will have non-oscillatory solution only if 


4 
— <0. 4.32 
The above relation can also be expressed in terms of particle densities as 
d dn'— dn& 
fA -_ Ee - =| | (4.33) 
do aL do do aL 


The electron density including the potential profile at any point is given by 


1 + er f e(¢—¢o) 


€ € ed k BTy 
n°(@) = nyexp : (4.34) 
kT 1 + er f me 
ae 
Differentiating the above equation with respect to @ yields 
e(¢—b0) kpT¢ 
dn°($) e née es ed 1+ erf kpTy nF e €Po et (d—¢o) (4 35) 
do kT kel 1+ erf mate kel kpTy 1+ erf se 
But at z — L, 6 > 0, Eq. (4.35) reduces to 
dné e That; 
nN 2) = ies 1 lex eee et (¢o) ; (4.36) 
do j,4, Kel KpTy | 14 er f iat? 
The ion density at any point x is given by 
i Yo 2eo = i 
aG= dis c a | # (0); (4.37) 
aa mv 
After differentiating the Eq. (4.38) we get 
dn'() / 8 FE) 2e¢ |? 
= du—— |1- —> | . 4.38 
dp oe ergs mv ee) 


As 2e¢ < m'v? at x & L, we can expand above equation to get 


fiw _ 3e¢ 
[| - a, dv a ae aa. (4.39) 
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AtzeL 


pepe 


From Eq. (4.33), (4.36) and (4.40) we get 


‘ edo 
(ay [eo 
Ud gf Pudge eT Do l+terf nat 


where 
ci ora [va air) 
UF ge Mie Is v2 
and 
Ce — kpT; 
s m % 


C’, is ion acoustic velocity. 


4.4 Presheath-sheath approximation 


32 


(4.40) 


(4.41) 


(4.42) 


(4.43) 


(4.44) 


The plasma flowing towards the wall passes through the two regions: the narrow region in 


which there is large gradient of electric field and supersonic velocity called “sheath”, an 


d 


the region attached with the bulk plasma with relatively weak gradient of variables and in 


general subsonic flow is called “presheath” [20]. 


The scale length of the sheath and presheath is different. So usually they are studied sepa- 


rately using different models and methods [6]. For our case, to couple the different solutions 


at the presheath and sheath side is important. 
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we have applied the coupling scheme [16], which gives the plasma parameters at the sheath 
side of the sheath entrance for parameters given at the presheath side. This coupling scheme 
satisfies the most crucial requirement of the presheath-sheath transition, i.e. quasi-neutrality, 


the sheath edge singularity condition and kinetic “Bohm-Chodura criterion”. 


4.4.1 Presheath 


We assume that the parameters at the presheath side of the sheath region are given or can be 


obtained by velocity distribution function. Let us assume the ion and electron densities, oe 


e e 
and n ps? 


. eye 4 
ps» the fluid velocities, u,, and u 


and the temperatures T’ Le and T° ,; at the presheath 
side of the sheath-presheath boundary (characterized by x = L). These six parameters are 


not completely independent but are related by the quasineutrality condition 


n,=n (4.45) 


ui, Bes (4.46) 
where 
k ‘Tt + it Ned 
ie uae bs) (4.47) 
m* 


Also the electric current density at presheath region is defined by 


Jes =e (n' ul — nN, ,U, ) (4.48) 


ps “ps ps "ps 


and at the presheath region 
Ns = Nis = ips (4.49) 


In summary, at the sheath edge (a = L) our two fluid presheath plasma is characterized by 


e 
ps? 


nine parameters, Nis; Nps: Nps: Ups: Ups: Ls, Lys, C® and jys, which are related by conditions 
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(4.45)-(4.49) so that knowing four of them we obtain the remaining five. 


4.4.2 Sheath 


The electron and ion distribution function at « = L describing the sheath are given by Eq. 
(4.8) and Eq. (4.10) respectively and first three moments are calculated from Eq. (4.14), Eq. 
(4.19) and Eq. (4.23), for electron and from Eq. (4.17), Eq. (4.21) and Eq. (4.25) for ion 
respectively. 

The seven sheath parameters indicated in sec. 4.2(a) are no independent of other but are 


interrelated by two boundary condition at x = L, namely quasineutriality condition 


ny (A®, Ts, Po) = ny (A*, Ts, Vines Vor) (4.50) 


and the Bohm Criterion, 
vee dy mi(v—viir)* nem —kT; De® 
—e = ae = -(1 4 4.51 
i pel oer = ered! + Gag, OP oe 


Hence, prescribing five of seven parameters will fix the remaining two, Eq. (4.50) and (4.51), 


following from the Bohm criterion. 


4.4.3. Presheath-Sheath coupling 


Now,we need to couple our kinetic sheath solution x = L_ to the fluid presheath solution 
at « = L,. For this, we require the sheath side particle density, nj, fluid velocity, uf, and 
effective temperature, T?;,, respectively to be equal their presheath side counterparts, nyps, 


Upe, BNA Te. 
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Equating these values taken from presheath counterparts, we obtain, 
TkTe 

Nis = ae C* A® 
kT Dé 
Und =A} — 
ps Tme Ce 

Aes. DP - 2. CPs 

T°. =Tell—,/- — 
ps = 7 \) «kT? Ce x De? } 


ti f ri yi 
Nps = |! ori C"A 
ib us 2kT} D 
ps “mL ami Ci 


and 
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(4.52) 


(4.53) 


(4.54) 


(4.55) 


(4.56) 


(4.57) 


Hence, we have seen equations, Eq. (4.51) and Eq. (4.52-4.57), which permit us to express 


seven sheath parameters A®, Tf, bo, A’, T;, Uj,z, and v},, in terms of five presheath parame- 


e e a a 
ters, Nps, Ups, Ips: Ups and Ty. 


Thus, for given presheath parameters the correspondinfg sheath parameters can be obtained 


by solving these equations. 


Chapter 5 


Numerical Method 
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5.1 Discretization of the simulation region 


Our simulation region x = 0 and x = J, is discretized uniformly in configuration and velocity 


space as shown in Fig. 5.1. 


Figure 5.1: Schematic representation of the phase-space griding. The vertical and horizontal 


dashed lines respectively represent the position and velocity grids. 


In Fig. 5.1 the position grid point is denoted as x; (j = 1,2,3....... ,Nx), with n, the total 
number of grid points, and the separation between consecutive grid points is denoted by Az. 
In our discretization, 7 = 1 and 7 = n, respectively correspond to the left-hand boundary 


(x = 0) and the right-hand boundary (a = L). In our simulation we choose n, such that the 
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grid size is always less than the Debye length of the injected electrons 


€ok ple 
a \| ; eae 


If we are interested to estimate the values of any quantity that we required in simulation 


region we use method of interpolation. Let the value of any quantity Q(x) at the grid points 
are denoted as Q; = Q(x;), its value at any point ‘x’ between two grid points Q; and Qj41 


can be approximated by means of linear interpolation approximation method 


ate) = (3) a+ (GB) ain (652) 


where2 7 P<. My4. 


Another standard approximation is the second-order interpolation between the points x j—1 


and x41 is 

Q(x) = Qj +a(x— 23) +(e — 2), (5.3) 
at point 7 = xj;_1 1s 

Qj-1 = Q; — aAz + b(Az)? (5.4) 
and at points x = 2 ;41 is 

Qja1 = Qj t aA + b(Az)?. (5:5) 


From Eq. (5.4) 


ao (4 a Git er) (5.6) 


Putting the value of ‘a’ in Eq. (5.5), we get 


Qj41 = Qj + [Q; ae eRe b(Az)?] (5.7) 


Solving above equation 


_ (G2 te t) 
2(Az)? 
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Using this value of ‘b’, we get value of ‘a’ as 


5 (Gu ee) . (5.9) 


Putting the value of ‘a’ and ‘b’ in Eq. (5.3) 


ay aetenay (Qh) Year ( WAM) 


The Taylor series expansion about 29 is 


(x — x)? 


SO" 0) ene (5.11) 


Q(x) = Q(ao) + (# — 20) Q' (x0) + 


Comparing Eq. (5.10) with Eq. (5.11), the first and second derivatives of any quantity Q(x) 


are approximated respectively as, 


dQ\ — (Qi4i — Qj-1 
(42) - (929) aa 


GEO Ve i OQ paa ae Oo, 
(aa),= ( a) ms 


5.2 Electron density distribution 


and 


In simulation region the electron density calculated analytically. The electron density at 
some grid point x; is obtained by using the electron distribution function in the expression 
for particle density Eq. (4.11). In order to obtain the electron density profile at any point x 
we require potential ¢(x) to be known. In our calculation, we obtain the electron density at 


the grid point x; in terms of potential as 


Bn ine ed; 
ni = nyexp aT 


where n§ = n°(x;) and $; = $(2;) . 


l+erf ee (5.14) 
f 


With the help of Eq. (5.14) electron density at any point in simulation region can be obtained 


if once we know the value of potential at that point. 
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5.3 Discretizing ion velocity space with a fixed grid 


As in the Fig. (5.1), the ion velocity space is discretized uniformly with n, the total number of 
velocity grid points and the k“" ion velocity grid value being denoted as vz (k = 1, 2,3, ...., Nv). 
In our discretization k = 1 and k = n, respectively corresponds to the fastest and the slowest 
ion velocities taken into account. Since in our simulation scheme only ions with negative 
velocities are considered, we set vp, = 0 and choose | v,; | sufficiently large such that the 
velocity of the fastest ion reaching the left-hand boundary does not exceed this value, i.e. 
In our case we choose n, to be large enough for the velocity grid size to 


| U1 |2| Or aad |. 


always be considerably less than the ion thermal velocity Ui f- 


5.4 Jon trajectories 


5.4.1 Discrete set of the injection velocities 


In our simulation region the ions are injected at the right-hand boundary i.e at x = L in 


the velocity range —oo < v < v',. In our numerical implementation, we approximate this 


a 


rie be OS vi, and discretize the velocity domain with nj-q equidis- 


domain by the finite v 
tant ion injection velocities v'; (7 = 1,2,...,Ntra), Where T = 1 and T = Nira respectively 
corresponds to the injection velocities vj,,,., and vi,. Choice of vi, depends on the problem 


a 


mar, 18 Chosen sufficiently large for the ion 


considered and the maximum injection velocity v 
distribution function corresponding to velocities larger than this velocity to be negligible. In 


a 


our computations, we use Ujax., = Ve, — 4U;¢, Where vj, is the ion thermal velocity. 


5.4.2 Discrete set of ion trajectories 


The ion injection velocity in the simulation region at right hand boundary (x = L) represents 


the starting velocity of a collisionless ion trajectory, the irq discrete injection velocities vu! 
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define n, related trajectories labeled by the index 7. We denote the trajectory corresponding 
to some injection velocity v',; by vi(Z), and the velocity with which it crosses a grid point 
a; is UL; = vi (a;), so that vi, & vt, = vi (ZL) in phase space. The velocities vi represents 
the intersections of ion trajectories v = v},,, (7) and spatial grid lines x = xj, which is said to 
be “intersection velocities”. The ion injection velocities v!,; are independent of the ion grid 
velocities vz, but in practice we will often choose them such that each v', coincide with one 


of the v;’s. 


5.4.3. Numerical calculation of ion trajectories 


In the KTS method we trace the collisionless trajectories to calculate the related velocity 
distribution functions along them. In this section we consider only ion trajectories and hence 
omit the species index ‘i’. For the calculation of the ion trajectory, we simply discretize the 


ion equation of motion Eq. (3.13) and Eq. (3.17) in a time-centered manner as 


gts <= i 
= y™. 5.15 
=u (5.15) 
For change in v; 
vm — ym 2 Pia: 
At ‘ 
é€ 1 m—t 
= — |E(a™ 1) + (vx B), | 
m* 
= (2-4) (5.16) 
m* 


Where we use At is the numerical time-step size and m > 0 is an integral or half integral 
superscript, i.e. v™ = v(t™), amtz =o Gas etc. This corresponds to the total time 
elapsed in steps of At for an ion having been injected at time zero at the right-hand boundary. 


So, m = 0 corresponds to t = 0 with the injection (or starting) values 
and 


eae (5.18) 
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The choice of positions at half integral times and velocity at integral times in Eq. (5.15) and 
(5.16) makes the numerical scheme time centered. For any time step m > 1, with x~? and 
v™—! being given from the (m — 1)** step, we calculate the new ion velocity from Eq. (5.16) 


as 
2) (5.19) 
and we can also find new position of the ions from Eq. (5.15) 

gett = g™2 + Atu™. (5.20) 


The electric field E(x) in equation (5.19) at any point in the simulation region is obtained by 


the relation (3.28). Thus the numerical scheme ends up with x-points at half-integral times 


m 


™\ along a collisionless trajectory. We can also 


(ars) and velocities at integral times (v 
find the ion velocity at gma by the relation 


m+4 = ha + oa 


4 ; (5.21) 


5.4.4 Intersection velocities 


The velocities at half integral time as in Eq. (5.21) may not coincide with our fixed grid 
points z;. For the calculation of the intersection velocities for the 7 trajectory at some 
inner grid point x; the linear interpolation can be used as 
al pen 
Caz — x3) Up e+ («; — ata) Uz * 


a 5.22 
TJ Cae. — att) ( ) 


where x”-2 and x™*? are the adjacent points along the trajectory on the right-hand side and 
the left-hand side of the grid point x,;, respectively. In order to calculate the ion velocity at 


the left-hand boundary (x=0) we use a parabolic extrapolation scheme. For this considering 
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the last three calculated points along the trajectory which lie to the right-hand side of the 


1 3 


left-hand boundary [eee Un ) : Ca a) and (oe%, Un bs | , the ion velocity near 
the left-hand boundary using the parabolic extrapolation scheme can be approximated as 
m—2® 5 5 2 
UL) =v, * +a(x—-a™#) +6(z-a™#) ; (5.23) 


where ‘a’ and ‘0’ are constants. The value of the constants are calculated. At the points 


C= 3, the velocity is approximated as 


m—24 m—2® 5 2 
Ur 7 =Vz 7 +4 (a3 = at) +b Cae — x) (5.24) 
Solving above equation, we get 


1 5 2 
m—>5 m— > 1 eee) 
Uz * — Uz 7b (amb — a” s) 


i ( =a =H) (5.25) 
x 2—f 2 
At the point x = r™-2 the velocity is approximated as 
m—32 m—2® 3 5 3 5\2 
Un SE Ae (a — a3) +b (a3 — x) (5.26) 


Using the value ‘a’ from Eq. (5.25) in Eq. (5.26), we obtain 


mit 


m—® 5\2 
5 Uz 7 — Ug ?—b(a” 2—_™ s) ( 
| 


5 _3 m—4 m2 _5 i m—3 m2 
Ca -_ ae ° ) (ur _ re ; ) = Ga oe a | (ur oa ve ; ) 
b= : : ; (5.28) 
2. 


ake 1 5 2 3 5 
22230. _3 m—s5 m—s 2/5) Lak m—s m—s 
(2” aac = ; ) & ; 7 a ’) 7 (« a a: ° ) (ur ° 7 te ° ) 
= . (5.29) 


5 1 5 3 3 it 


Then using the value of constants ‘a’ and ‘b’ from Eq. (5.28) and (5.29) onto the Eq. (5.23), 


the ion velocity at the left-hand boundary is calculated as 


5 5 2 
?—arz™"2+b (2-3) , (5.30) 
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5.5 Ion velocity distribution function 


When we trace the ion trajectories for all nq injection velocities v',, we can calculate the 
corresponding intersection velocities Eq. (5.22) and the related distribution functions at all 
grid points x;. The injection velocity from where we start are usually uniformly spaced, but 
the corresponding ion velocities at the other grid points, 0 < x; < L, are non-uniformly 
spaced and usually do not coincide with the fixed grid velocities, v,. For our model, it is not 
necessary to know the ion distribution function at these uniformly spaced grid velocities. 

We can calculate the ion distribution function at the fixed grid points (position and velocity) 
by linear interpolation of the values associated with the trajectory intersection points to the 


nearest desired grid points [16]. 


5.6 Ion density distribution 


The ion density distribution n‘(x) is obtained from Distribution Function Approach (DFA). 
In the DFA, the ion density n'(x;) is obtained by integrating the ion velocity distribution 
function over velocity space. For the plasma model, we obtain the ion distribution function 
only at the intersection points [x;,v'(x;)] and can also be obtained at the fixed grid points 
[x;, Uz] where the distribution function is given. The ion density at x;, where the distribution 
) and the intersection velocities of all nq ion trajectories vs, are given, 


functions f'(x;, v1, 


may be obtained by the following discretized form: 


. ae . ; 
n'(xj) = 5 a [F'(x5, 05) + fey 4a g)) (Oraag — Urs) (5.31) 
T=1 
ieee esc, 1% or ; 
n'(L) = 2 [f°(L, v4) za f(L, vr 41,1) (Caen > vir) ; (5.32) 
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5.7 Solution of Poisson’s equation 


After the calculation of the electron and ion densities respectively from Eq. (5.14) and Eq. 


(5.32), the space charge density distribution can be obtained from 
p(z3) = S° g’n*(x;) (5.33) 
= e[n'(xj) —n°(x,)]. (5.34) 
By the knowledge of the space charge distribution p(x,;) from Eq. (5.33), we can solve the 


Poisson’s Eq. (3.29) numerically to obtain the related new potential distribution, ¢(x;). 


Using Eq. (5.13), the discretized form of Poisson’s equation is given by 


bj41 — 20; + O5-1 Pj 
=——, 5.35 
(Ag? Eo oP) 
Now we write this equation for all internal grid points in the simulation region (7 = 2,3,....,N2— 
1) which yields the following n, — 2 equations involving the n, unknowns, @}, ¢,...., On,: 
Az)? 
bi — 22 + 3 = uk ; po (5.36) 
hae 
b2 — 263+ b4 = Erp, (5.37) 
_ (Az)? 


a (5.38) 


Pn —2 —_ 2Prig—1 oF One Z 


We have fixed the potential values at the two boundaries as 


o(t=L) = on, 
= OL (5.39) 
and 
d(x =0) = G1 


= do. (5.40) 
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By solving the Eq. (5.36), Eq. (5.39) and Eq. (5.40) we now obtain the potential distribution, 


i.e. the potential values ¢1, ¢2,....,@n, Which can be expressed as a single matrix equation as 
TO! 50520 0 oi o 
i -203- ti 60 0 oe — ON py 
0 1 -2 0 
= (5.41) 
r 2 
1-21] ] bng-1 in, 
O08. One OL 


The matrix on the left-hand side and the one at the right-hand side are known. In order to 
solve the matrix equation inverse of the first matrix is multiplied with the right-hand side 


matrix which is done using a simple command in MATLAB program. 


5.8 Relaxation scheme 


The exact solution of Poisson’s equation converges only for short systems (a few Debye 
lengths). As the system length increases, small fluctuation of the potential causes unphysical 
accumulation of the charges and the scheme breaks down. To overcome this difficulty we 
use the relaxation scheme [21]. The numerically exact new potential distribution om (xj) is 
obtained by numerically solving Eq. (5.41) which is linearly combined with the old potential 
distribution ¢°"—) (x;) to obtain the “re-adjusted” new potential distribution ¢")(2;), which 


is actually used as the relevant result of the m*” iteration: 
gp (x5) = wo (a) + (L- wd (a), (5.42) 


where w is relaxation parameter, having with 0 < w < 2. 
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5.9 Electric field calculation 


Here, we will discuss how the electric field is calculated at all the grid points from a given 
potential profile, and how the field is calculated for the points lying between two grid points. 


(a) Electric field at the inner grid points: 


At the inner grid points j = 1, 2,.......... , Nz — 1, the electric field is obtained as 
(x 541) — O(2j-1) 
E(x) =—- 5.43 
(7) 2Az ( ) 


(b) Electric field at the left-hand boundary: 
Near the left-hand boundary (x = 0) we assume, for the potential, the parabolic approxima- 


tion as 
d(x) = ¢1 +x + ba’, (5.44) 


where ¢; = @p is the potential at x = 0. 


For the second and third grid points, we have the two equations 
bo = d1 + aAz + b(Azx)? (5.45) 
and 
63 = $1 + 2aAz + 4b(Az)? (5.46) 


respectively, from which the constants ‘a’ and ‘b’ can be obtained as 


_ 361 — 462 + 3 
=— AG (5.47) 
and 
= _ 1 = 22 + 3 (5.48) 


2 Aw)? 
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The electric field near x = 0 is obtained by taking the derivative of Eq. (5.44) as 


E(x) = a (5.49) 
eine (5.50) 


So that the electric field at the left-hand boundary is given by 


E(x =0) = -a (5.51) 


301 — 462 + 3 


52 
2Az 40.82) 


(c) Electric field at the right-hand boundary: 
Near the right-hand boundary (x = L), we assume for the potential the parabolic approxi- 


mation as 
(x) = bn,—2 + C(# — &n,-2) + d(a — tn,-2)’, (5.53) 


where ¢y,—2 is the potential at the grid point r,,-2 = L — 2Az. For the grid points ry,-1 


and %,,, we have the equations 

Pne—1 = Png—2 + CAx + d(Az)? (5.54) 
and 

One = Ong—2 + 2cAzx + 4d(Azx)” (5.55) 


respectively, from which the constants ‘c’ and ‘d’ can be obtained as 


Oa = 4@n,-1 + 30H 20 


pa mG (5.56) 
and 
On — 2¢n -1 a On —2 
= “ “ Ot 
4 2( Ae)? One 
The electric field near x = L is obtained by taking the derivative of Eq. (5.53) as 
do(x) 
E = 
(x) a 


= —c—2d(x—%p,_,), (5.58) 
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so that the electric field at the right-hand boundary is given by 


—c — 4dAx 


— 80nx = 4dng—1 + Oita3 
2Ar 


SS 
5 
| 
oS 
| 


(5.59) 


(d) The electric field at any point in the simulation region: 


Once we have calculated the electric field at every grid point, the electric field at any point 
in the simulation region is obtained by linear interpolation between two nearest grid points 


x; and Xj41 as 


(2j41 — ©) B; + (@ — 25) F541 


(zs) = eG 


(5.60) 


5.10 Iteration scheme 


We are dealing with collisionlal, time independent system only. The initial-guess potential 
$(x;) is taken as input to the main Iteration Block (located between the points P and Q 
in Fig. (5.2), which will yield the first Iteration to the potential distribution 6 (x,;). With 
the new potential as input, the main Iteration Block will be invoked again yielding ¢°)(2;) 
until the potential distribution has converged in the sense outlined in Sec. (5.11) below. 


(a) Boundary conditions: 


The boundary potentials ¢(z = 0) and ¢(a = L) and the boundary injection distribution 
functions f*(L,v) are to be provided as input parameters. Hence, they must be specified 
before entering the iteration scheme and are kept constant throughout the entire simulation. 
(b) Initial guess to ¢(z): 

To start the scheme, we must suitably prescribe an initial-guess potential distribution. As we 


restrict ourselves to potential distributions ¢(x), which decrease monotonically from ¢(L) = 0 
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to ¢(0) = @ < 0. The starting potential distribution is chosen to be a linear interpolation 


between the potential values at the boundaries. 


5.11 Main iteration block 


Here we describe the main iteration scheme discussed in the Fig. (5.2). The main iteration 
block carries out the m“” iteration [i.e. it calculates the new distributions ¢”(zx;) | for a given 
input (old) potential distribution ¢"~'(x;) by performing the following three steps: 

Step 1: The new electron density distribution n° (x;) is calculated analytically using 
Eq. (5.14). The new ion density n‘(zx;) is calculated by the distribution function ap- 
proach(DFA). In DFA, the new ion densities are obtained by velocity-space integration of 
the new velocity distribution functions |cf. sec. 4.6 (a)]. 

Step 2: From the new densities n*(™ (x;), obtained in step 1, the new space-charge density 
pe’ (x;)is calculated using the Eq. (5.33). 

Step 3: The new potential distribution ¢°”)(x;) is calculated by solving the matrix Eq. 
(5.41) numerically, with the new space-charge density p” (x;) inserted on the right-hand 


side. 
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Specify boundary potentials ( @ and @, ) and injection distribution 
functions f;(Z.v)- 


Prescribe initial-guess potential distribution ¢°” ( x; ) 


Step 1: 
Calculate n* (x ) by means of the analytical approach. 


J 
Calculate 1‘ (x ) by the distribution function approach 


Step 2: calculate p™”’ (x, ) 


Step 3: find ¢”’ (x,) by solvingthe matrix equation. 


Check of convergence 


lo" (x )-9 (x) <6 


yes 


Output, the self-consistent quantities 


9° (x)= 0 (x) 0° (x).0 (x), 


Figure 5.2: The iteration scheme. 
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5.12 Convergence check 


The new potential distribution, 6°” (x;), obtained in Step 3 of the main iteration block is 
compared with the old potential distribution, ¢°"~) (x;), and we consider the convergence 


to be reached if at each point x; the condition 
| o™(x;) — 9 Y(a;) |< 6 (5.61) 


is satisfied, where 6 is defined as accuracy parameter. 


5.13 Numerical parameters 


The Kinetic Trajectory Simulation is applied to our problem. We specifically consider the 


following parameters as given at the presheath side of the sheath-presheath boundary: 
e Hydrogen plasma with plasma density (n,,) = 10'8 m~*. 
e Electron temperature (T%,) = 10eV, ie. 10 x 11604.9 K. 
e Ion temperature (T/,) = 0.1T¥,eV. 
e Ion polytropic constant (7*) = 3. 


The resulting dimensional sheath parameters at the sheath side of the sheath-presheath 
boundary are then obtained from the coupling scheme. We choose the sheath to be of width 
10 A%, where Af is the electron Debye length at the sheath entrance associated with the 
injected electrons. 

The discrete ion injection velocities at the sheath entrance is discretized uniformly with 300 
injection-velocity grid point such that the ion injection velocity grid step is considerably less 
than the ion thermal velocity. The region between x = 0 to L is discretized uniformly with 
31 grid points. We define the system to have reached convergence (cf. sec. 5.12) if the 


maximum difference in potential values before and after iteration equals 10~° V or less. 


Chapter 6 


Results and Discussion 
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In this section, the results obtained by simulation of 1dlv plasma sheath by introducing 
collision term which is characterized by a [cf. Section 4.1], are presented and discussed. The 
value of a represents the degree of collision, i.e. how frequent the collision occur. Lets take 
a = 0.01 which implies that if a particle encounters 100 other particles it undergoes a single 
collision, a = 1 corresponds to the case where every encounter is collision . For collisionless 
case a = 0and 0 < a < 1 for collisional case. The plasma sheath structure have been studied 
for different values of a, for the same presheath parameters. The final self consistent profile 
for ion and electron densities, total charge density, potential distribution, average velocity, 


and effective temperature within the sheath region have been plotted and discussed. 


6.1 lon density profile 


Figure 6.1: Self consistent ion density versus normalized distance for different values of a 


The Fig. 6.1 shows the self consistent ion density versus normalized distance for different 
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values of a (degree of collision). In this chapter in all plots the distance from the wall 
is normalized with respect to electron Debye length at the sheath entrance and is termed 
normalized distance (§ = #/Ap). It is observed that the ion density decreases monotonically 
from the sheath entrance towards the wall acquiring its minimum value. The decrease in ion 
density is due to electrons have higher thermal velocities, reach the wall faster and hence 
the wall will be at negative potential. This repels further electrons thereby lowering their 
densities in the sheath region. This in turn causes the ion density to decrease conserving the 
flux. The plot for a = 0, is shown to compare our result to previous work in collisionless 
cases. At sheath entrance the ion density is maximum and minimum at wall. As the degree 
of collision increases the rate by which ion density changes is different, starting with the same 
value of presheath ion density. At any given point in the sheath region the ion density is 
larger for higher value of a, the ion density is minimum for collisionless case, i.e., a = 0, 
increases as we increase the value of a which explained on the fact that collision reduces the 
velocity of ion for flux conservation it’s density must increases. 

Increasing the value of a the rate at which the ion density increase is decreased and after 
certain value of a the ion density becomes saturated. Increasing the value of a from 0.2 to 
0.3 the rate of increase of ion density is not sharp as compare qa increases from 0 to 0.005. 
As we move towards the wall the ion density started decreasing at distance approximately 
10\p from wall for a = 0 (collisionless case), between 10\p and 9Xp for a = 0.0025, between 
9X\p and 8Ap for a = 0.01, between 8Ap and 7Ap for a = 0.02, between 7Ap and 6Ap for 
a = 0.03. Increasing the value of a the point after which the ion density starts to decrease 
rapidly shift towards the wall. The similar ion density profiles can be obtained for higher 


value of a but the profiles doesn’t changes significantly. 


6.2 Electron density profile 


The Fig. 6.2 shows the self consistent electron density versus normalized distance for different 


values of a. The electron density decreases monotonically from sheath entrance to wall. Since 
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Figure 6.2: Self consistent electron density versus normalized distance for different values of 


a 


electron density obtained analytically which explicit dependent on potential, as we fixed the 
wall potential constant causes the electron density constant what ever be the value of a. As 
we increase the value of a the electron density increases at any point in sheath region for 
flux conservation due to increase of ion density. At sheath entrance the electron density is 
maximum and minimum at wall. As the degree of collision increases the rate by which electron 
density changes is different, starting with the same value of presheath electron density. As 
increasing the value of a the rate at which the electron density increase is decreased and 
after certain value of a the electron density becomes saturated. Increasing the value of a 
from 0.005 to 0.3 the rate of increase of electron density is not sharp as compare q@ increases 
from 0 to 0.005, the electron density increases sharply. Increasing the value of a the point 
after which the electron density starts to decrease rapidly shift towards the wall. The similar 
electron density profiles can be obtained for higher value of a but the profiles doesn’t changes 


significantly. 
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6.3 Comparison between ion and electron density pro- 


file 


Particle density (m73) 


Figure 6.3: Self consistent ion and electron density versus normalized distance from the wall 


for different value of a 


The Fig. 6.3 shows self consistent ion and electron density profiles both for collisionless 
(a = 0) and typical collisional (a = 0.3) case. In both case ion density is greater than electron 
density in sheath region and at wall but equal at sheath entrance (to satisfy quasineutrality 
condition). The electron density decreases more rapidly than ion density. Due to the negative 
potential of the wall, the electrons are reflected and hence the density of electrons as well as 
that of ions density decreases to conserve flux also the ion thermal velocity is less than that 
of electrons, the ion density exceeds the electron density so that the flux is conserved. For 
better comparison we can compare the density profile for different value of a as plotted in 


Fig. 6.1 and Fig. 6.2. instead of taking typical value of a. 
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6.4 Potential profile 


Figure 6.4: Self consistent potential versus normalized distance from the wall for different 


value of a 


The Fig. 6.4 shows self consistent potential profile for both collisionless (a = 0) and typical 
collisional (a = 0.03) case. In both case negative potential increases monotonically from 
sheath entrance to wall. As we fixed the potential at wall and at sheath entrance fixed in 
simulation, the potential does not change at wall and sheath entrance but in other sheath 
region the potential has high value for high value of a. Also increasing the value of a the 
negative potential starts to increase at the point nearer to the wall. The similar potential 
profiles can be obtained for higher value of a but the profiles doesn’t changes significantly. 


For more exact potential profile it is better to plot graph for different value of a. 
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6.5 Electric field profile 


Figure 6.5: Self consistent electric field versus normalized distance from the wall for different 


value of a 


The Fig. 6.5 shows the self consistent electric field profile versus the normalized distance 
from the wall for both collisionless (a = 0) and collisional (a = 0.03,0.005) case. The 
electric field is always negative in the sheath region due to the negative potential and its 
magnitude increases towards the wall. Electric field is maximum at wall and minimum at 
sheath entrance. The field at wall is higher for higher value of a and becomes saturated 
as we increase the value of a due to high ion density but same electron density. Increasing 
the value of a from 0.005 to 0.3 the rate of increase of electric field is not sharp as compare 
q@ increases from 0 to 0.005, the electric field increases sharply. Increasing the value of a 
the point after which the electric field starts to increase rapidly shift towards the wall. The 
similar electric field profiles can be obtained for higher value of a but the profiles doesn’t 


changes significantly. 
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6.6 ‘Total charge density profile 


Figure 6.6: Self consistent total charge density versus normalized distance from the wall for 


different value of a 


The Fig. 6.6 shows the typical self consistent charge density distribution versus the nor- 
malized distance from the wall for different value of a. The total charge density is zero at 
the sheath entrance (quasineutrality) and increases as we move towards the wall, where it 
acquires it’s maximum. As the value of a increases the charge density at wall increases due 
to increase of electron and ion density. Increasing the value of a the point from which charge 
density rapidly start to increase shift to the wall. The similar charge density profiles can be 


obtained for higher value of a but the profiles doesn’t changes significantly. 
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6.7 Average ion velocity profile 


The Fig. 6.7 shows the self consistent average ion velocity versus the normalized distance 
from the wall for different value of a. The velocity of ion increases as we move from sheath 
entrance to wall in both collisional and collisinaless cases. As the value of a increases the 
velocity of ion decreases due to collisional effect. The collision always decrease the velocity 
of colliding particle that same happen in the sheath region. The similar velocity profiles can 


be obtained for higher value of a but the profiles doesn’t changes significantly. 


Figure 6.7: Self consistent ion average velocity versus normalized distance from the wall for 


different value of a 
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6.8 lon effective temperature profile 
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Figure 6.8: Self consistent ion effective temperature versus normalized distance from the wall 


for different value of a 


The Fig. Fig. 6.8 shows the self consistent ion effective temperature profile versus normalized 
distance. The effective temperature decreases from sheath entrance to wall which is due to 
decrease in ion density towards the wall. As collision increases the effective temperature at 
any point in sheath region also increased due to increase in ion density by collisional effect. 
The similar temperature profiles can be obtained for higher value of a but the profiles doesn’t 


changes significantly. 


CHAPTER 6. RESULTS AND DISCUSSION 63 


6.9 Discussion and conclusion 


We have used a 1dlv kinetic model for the study of plasma sheath structure for different 
values of a (degree of collision) keeping other parameters the same. In all cases, the physical 
parameters such as potential, ion density, electron density, total charge density etc have very 
weak gradient close to the sheath entrance, however, they all have sharp gradient close to 
the wall. 

The ion and electron density increases at any point in sheath as the value of a is increased. 
The ions are then accelerated towards the wall and are absorbed by the wall, decreasing the 
negative wall potential. This causes the increase in the electron density. However the system 
balances itself for higher n; than n, to conserve the flux (electron have higher mobility than 
ion) for all valid values of a. 

Both the ion and electron density decrease gradually from the sheath entrance towards the 
wall. Precisely, the density of electron decreases faster than that of ion. This is due to the 
reflection of electrons by negative potential at the wall. The decrease in ion and electron 
density is less prominent as we increase the value of a at the sheath. 

All sheath parameter measured above have great effect of a. If we compare the effect of a 
when its value close to 0 and for higher value, it is observed that when value of a are varied 
in lower region close to 0 effect the sheath structure highly than varying the value of a in 
higher region. From our simulation result we can say after certain value of a, the collision 
does not significantly effects the sheath structure. 

The collision effect decrease the sheath width which can be explained by potential profile. 
The sheath entrance starts from where negative potential starts increase. As degree of colli- 
sion increases the negative potential starts increase at the point nearer to the wall. This also 
explained by the total charge density. The plasma region is quasineutral region but sheath 
region is not quasineutral. From which the quasineutrality breaks we can say the sheath 
started. 


The results discussed above suggest that sheath structure is greatly influenced by collision 
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at the sheath region. The work developed here starts a suitable basis for an important study 
in future, the main purpose being to arrive at more realistic collision models. Besides, our 
work is expected to provide necessary idea for the implantation of ions at material walls as 
well as in various applications of plasma. 


The following cases may be considered in future work: 


(a) Time dependent cases 

(b) Magnetic field into consideration 

(c) Higher dimensional analysis 

(d) Considering the specific collision 

(e) Considering time dependent collision term 
(f) 


f) Kinetic treatment to electron 


Appendix A 


Matlab Files 


In this section we present the MATLAB files. 


(a) input file 

alpha= input(‘Type the value of alpha’) ; 

% Basic parameters: 

kB = 1.38062e — 23; % Boltzmann constant 

epsilon = 8.85419e — 12; % permittivity of the medium 
e = 1.602192e — 19; % electronic charge 

Ze =e; % ion charge 

Mi = 1.672e — 27; % ion mass 

Me = 9.109e — 31; % electron mass 

mu = Me/Mi; % electron to ion mass ratio 

gammai = 3; polytropic constant for ions(isothermal case) 


% presheath parameters: 
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Tpse = 10 * 11604.873; % electron temperature at presheath side 

Tpsi = 0.1 * Tpse; % ion temperature at presheath side 

nps = 1e18; % plasma density at presheath side 

cpsi = sqrt(kg * (Tpse + gammai * Tpsi)/Mi); % ion acoustic velocity at presheath side 
upst = —cpst; 

% Normalized quantities: 

TpsiN = Tpsi/Tpse; % normalized ion temperature at preseath side 
upsiN = upsi/sqrt(2 * kp * Tpse/Mi); % normalized ion acoustic velocity 
% sheath parameters: 

JpsN =0; %Normalized current density 

phiL = 0; % potential at x = L 

rhoL = 0; % charge density at x = L 

L=‘10*« DLe’; % system length 

% numerical input parameters: 

ntra = 200; 

phi0x = ‘phi0 « (1 — X/L); 

nx = 31; % number of x-grids 

dt = le — 11; 

deltaphi = le — 5; % desired accuracy in the iteration 


w = .08; % signifies the relaxation parameters 
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nv = 300; % number of velocity grid 
phif Ni = linspace(—5, —0.001, 401); % initial guess (always negative) 


TdUmin = —5.8877254338; % lower limit of variable ‘Tau’. useless to keep below: -5.8877254338, 


since f — 0. 
TcLiNi = linspace(—5.88, 0,90); % initial guess 


asymTcLiN = —4.5; % limit below which we use the asymptotic expression 


(b) Distribution function for ions at injection 


% specify the injection ion velocity distribution 


function for specified velocity vx function df = Dfun(v,) global niL vmLi vtfi kB Tfi Mi 


vmaxLi veLi 


df == niL * exp(—((vx — vmLi)/vt fi).*)/sqrt(.5 * pi * kB * T fi/Mi)/(erf((—umarLi + 
umLi)/ut fi) + er f((—umLi + vcLi)/vt fi)); 


(c) Poisson solver 


% solves the poisson equation for given charge density and potential at the boundaries fixed 


to some constant value 

MR(1) = phi(1); MR(nz) = phi(nx); MR(2: nz — 1) = —dx? * rho(2 : nz — 1)/epsilon; 
phin = MLM R’; phin = phin; 

phin = w * phin + (1 — w) * phi; phin(1) = phi(1); phin(nzx) = phi(nz); 


(d) Electric field calculation 


Here we calculate the electric field to the given potential function 
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field = DF AE field(z) 


global Ef dx X nx 


i —— 


field = Ef (1); 


elseif x == X(nz) 


field = Ef (na); 


else 


f= ceuwsar): 


field =(EfG)*(XG+)—#) + EfG +1) * (a — X(9)))/da; 


end 


(e) Coupling of the sheath with the presheath 


% 


% 


% 


% 


% 


% 


% 


% 


% 


Coupling sheath with presheath. 

For given presheath plasma parameters (ups, Tpsi, Tpse, Jps, gammai, mu, etc.) 
this program 

calculates sheath plasma parameters (Ai, Ae, vcLi, vmLi, Tfe, Tfi, phi0, etc.) 
required for our 1dlv plasma-sheath simulation 


This proceeds as following: 


1) Provide the name of your input file, where all required input parameters are 
specified. 
2) Solves the electron irreducible equation to obtain ‘phifN’, and then the other 
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% 


% 


% 


% 


% 


% 


% 


% 


% 


% 


% 


% 


% 


corresponding 

electron parameters (TfeN, AeN, phiON, etc.) are calculated. 

3) Solves the ion irreducible equation to obtain “TaucLiN’, and then the other 
corresponding 

ion parameters (TfiN, AiN, vmLiN, vcLiN, etc.) are calculated. 

4) Once all normalized parameters are known, the dimensionless physical 


parameters are calculated 


using their respective normalizing equations. 

close all, 

clear all 

solving the electron irreducible equation for ‘phifN’ — — —start — —— : 
Calculating the lhs (Le) of the equation (independent of ‘JpsN’ and “TpsiN’): 


Calculating other parameters: 


for rphifN =1: length(phif Ni) 


phifN = phif Ni(rphif N); 


Ce=1+erf(sqrt(—phif N)); 


De = exp(phif N); 


TfeN = 1/abs(1 — sqrt(—4 « phif N/pi) * De/Ce — 2/pi * (De/Ce)?); 


Le(rphif N) = De/Ce * sqrt(T feN/pi); 


% 


lhs of the electron irreducible equation 
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end % for r=1: length(phif Ni) 
% Calculating the rhs (Re) of the equation (function of ‘JpsN’, ‘TpsiN’, ‘gammai’ and ‘mu’): 
Re = sqrt(mu) « (JpsN + sqrt(.5 * (1+ gammai * TpsiN))); % rhs of the equation 


phif N_sol = interp1(Le, phif Ni, Re); % Interpolating for the solution (phifN) for which 
Dbe= Re, 


% Thus obtained ‘phif N_sol' is the required value. 

if isnan(phif N_sol) 

disp(‘Could not find the solution. Change the input-range for phifN and start again.’) 
pause 

else 

phif N = phif N_sol 

Ce =1-+erf(sqrt(—phif N)); 

De = exp(phif N); 

TfeN = 1/abs(1 — sqrt(—4 * phif N/pi) * De/Ce — 2/pi * (De/Ce)”); 

AeN = 1/Ce/sqrt(T feN); 


% solving the ion irreducible equation for ’*phifN’— — —end — —— 


% solving the ion irreducible equation for ‘TcLiN’— — —start — —— 


% Calculating the rhs (Ri) of the equation (function of ‘phifN’, independent of ‘JpsN’ and 
‘TpsiN’): 


Ri = (sqrt(pi) + De/Ce/sqrt(—phif N))/TfeN; % Calculating the lhs (Li) of the equation: 
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rr =1; for r= 1s léngth(feLiN1) 
TcLiN =TcLiNi(r); 


if TcLiN < asymTcLiN]% using asymptotic expansion: er f(x) = 1—exp(—x?)/sqrt(pi)/x* 
(1 — 1/2/22 + 3/22/24 — 1.3.5/23/x6 +...) 


Ci = exp(—TcLiN?’) /abs(TcLiN)/sqrt(pi)*(1—.5/TcLiN?+.75/TcLiN4—15/8/TcLiN®+ 


105/16/TcLiN® — ...945/32/TcLiN™) + 10395/64/TcLiN? — 135135/128/TcLiN™ + 
2027025 /256/TcLiN 6 —34459425 /512/TcLiN *+...654729075/1024/TcLiN —13749310575 /204 
316234143225 /4096 /TcLiN* — 7905853580625 /8192/TcLiN?): 


else % if TcLiN < asymTcLiN 

Ci=1+erf(TcLiN); % if TcLiN < asymTcLiN 

end 

Di = exp(—TcLiN?); 

TfiN = TpsiN/abs(1 —2* TcLiN « Di/Ci/sqrt(pi) — 2 * (Di/Ci)?/pi); 
umLiN = —sqrt(0.5 * (1 + gammai « TpsiN)) + Di/Cix sqrt(T fiN/pi); 
a=vumLiN/sqrt(T fiN); % abbreviated for simplicity 


if TcLiN < —a % checking the integrability 


% Integrating: 

Tau = linspace(Taumin, TcLiN, ntra); % Discretizing for integration 
dTau = Tau(2) — Tau(1); % Width of the Tau-grid 

for au = 1 nire 


FnTau(rTau) = exp(—Tau(rTau))/(Tau(rTau) + a)?; 
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end 
LHS = dTau* (sum(F'nTau) — 0.5 * (FnTau(1) + FnTau(ntra))); 
Li(rr) = LHS/Ci/T fin; 
TcliN M(rr) =TcLiN; % storing the values for which the equation is integrable. 
rr=rrt+l,; 
end % if TcLiN < —a 
end 


for r= 1: length(phif Ni) 


% obtaining the solution by locating the point of intersection of ‘Li’ and ‘Ri’: 
Yplot(TcLiN M, Lt) 

Ahold 

% plot({min(TcLiNM) max(TcLiNM)],[Ri Ril],‘k — —’) 

TcliN,ol = interp|(Li, TcLiN M, Ri); % required solution 


% This method of obtaining the solution can be RISKY. The point where we have found 
‘Li = Ri’ is 


% the point of marginal validity of the Bohm criterion, which, in our case, is Li < Ri. 


% When we consider the marginal point it may be possible that because of numerical 


limitations it 
% lies in fact in the region, where the Bohm criterion is not satisfied. Hence, in order to 


% be sure we take a point very close to the zero-point but still lying in the negative side of 
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% sure side of (Li — Ri), so that we may say that Li < Ri, in place of Li < Ri. 

% We have tested it many times and concluded that it has also negligible effect to the 
% other ion parameters derived therefrom. a = find(Li — Ri < 0); 

% points where Li < Ri a = a(length(a)); 

% points where Li < Ri, a = a(length(a)); 

% points closest to the zero-point but still Li < Ri TcLiN_sol = TcLiN M(a); 

% value of the closest point to the zero-point — Li < Ri 


% plot (TcLiNM(a), Li(a),‘mo’) % solving the ion irreducible equation for ‘TcLiN’ — — 


end — — 
% calculating other dimensionless ion parameters: 
% Thus obtained 'TcLiN_sol' is the required value. if isnan(TcLiN_sol) 
disp(‘Could not find the solution. Change the input-range for TcLiN and start again.’) 
pause 
else 
TcliN =TcliN_sol 
% Calculating other parameters: 
if TcLiN < asymTcliN 
% using asymptotic expansion: 


er f(x) = 1— exp(—x?)/sqrt(pi)/ax * (1 — 1/2/x? + 3/2?/x* — 1.3.5/23/x® + ...) 
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Ci = exp(—TcLiN?) /abs(TcLiN)/sqrt(pi)*(1—.5/TcLiN?+.75/TcLiN4—15/8/TcLiN®+ 
105/16/TcLiN® — ...945/32/TceLiN™) + 10395 /64/TcLiN” — 135135/128/TcLiN™ + 
2027025 /256/TcLiN!°—34459425/512/TcLiN'8+...654729075/1024/T cLiN”°—13749310575/204 


316234143225 /4096 /TcLiN* — 7905853580625 /8192/TcLiN?°): 


else % if TcLiN < asymTcLiN 

Ci=1+erf(TcLiN); 

end 

% if TcLiN < asymTcliN 

Di = exp(—TcLiN?); 

TfiN = TpsiN/abs(1 —2* TcLiN « Di/Ci/sqrt(pi) — 2 * (Di/Ci)?/pi); 
umLiN = —sqrt(0.5 * (1 + gammai « TpsiN)) + Di/Cix sqrt(T fiN/pi); 
AiN = 1/Ci/sqrt(T fiN); 

ucLiN =umLiN +TcliN x sqrt(T fiN); 

end % if isnan(T'cLiN_sol) 

end 

% calculating dimensional physical parameters: 

% electron dimensional parameters: 

Tfe=TfeN x Tpse; 

utfe = sqrt(2*kB xT fe/Me); 

Ae = AeN *nps * sqrt(2 * Me/pi/kB/Tpse); 


neL_ = Aexvtfe x sqrt(pt) /2; 
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phi0 = phifN «kB xT fe/e; 

% ion dimensional parameters: 

Tfi=TfiN « Tpse; 

ut fi = sqrt(2* kB * T fi/Mi); 

Ai = AiN *nps * sqrt(2 * Mi/pi/kB/Tpse); 

niL = nps; 

umLi = umLiN * sqrt(2 * kB * Tpse/M1); 

ucLi = ucLiN * sqrt(2* kB * Tpse/M1); 

(f) Ion density in DFA 

Q = Zex dt/Mi; % defined for simplicity 

forT! =1:ntra % selecting an ion trajectory 

m=1; % initial time count 

VaL=VIi(T); % ion velocity for the selected trajectory at injection 
Vm(m) = val; % Starting velocity for ions at the injection(injection velocity) 


Xm2(m) = L+ .5*vrl «dt +Q * dt * DF AE field(L)/8; Ymoving by dt/2 for time cen- 


tering(new starting position) 
if Xm2(1) <= X(nz — 1); 
disp(‘Time step is too large. Press any key to continue.’) 
pause 


end 
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whileXm2(m) > 0andXm2(m) < L % solving the equation of motion until the ion crosses 


x=O0orx=L plane 
m=m-+1; % next time step 
Vm(m) = Vm(m — 1) + Q * DFAE field(Xm2(m — 1)); 
Xm2(m) = Xm2(m — 1) + dt * Vm(m); 
end 
forr=1:m-1, 
Vm2(r) = 0.5 * (Vam(r) + Vam(r + 1)); 
end 
At = [L,Xm2(1:m—1)]; Vt = [Vm2]; 


a= ((Xt(m—2)—Xt(m))?*(Vt(m—1) —Vt(m—2)) — (Xt(m—2) —Xt(m—1))?*(Vt(m) — 
Vt(m — 2)))/(Xt(m — 2) — Xt(m — 1))/(Xt(m — 2) — Xt(m))/(Xt(m — 1) — Xt(m)); 


b = ((Xt(m— 2) —Xt(m)) * (Vt(m— 1) —Vi(m — 2)) — (Xt(m — 2) — Xt(m—1)) *(Vt(m) — 
Vt(m — 2))) /(Xt(mn — 2) — Xt(m — 1))/(Xt(m — 2) — Xt(mn))/(Xt(m) — Xt(n —1)); 


vx0 = Vt(m — 2) +a* Xt(m — 2) +b* Xt(m — 2)": 


Xt = [Xt, 0]; Vt = [Vt, v0]; % Position (from x=L to x=0) and the corresponding velocities 


along the trajectory 
Vti = interp1(Xt, Vt, X); % Interploating for the velocity at fixed grid(x) 
Vil z= Vin: 
Is(T) = alpha * cpsi * L/(DLe)? * abs(V Li(T)); % collision integral 


Df (T) = dfun(VLi(T)) + Is(T); 
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clearVmV m2VtV tiX m2X tabmrvx0 


end 


% ni calculating with interpolating onto the fixed velocity grid points 


ni = zeros(1,nzx); 


fie ea Bee fm 


forT =1:ntra—-1 


ni(j) = ni(j) + .5* (DF(T) + Df (T + 1)) * abs(V,7(T + 1) — Vz2(T)); 


end 

clear V,iT 

end 

%% calculates average velocity for ions at any point X 
ui = zeros(1,nz); 

JOR) Slee: 

Vit = Vit, 9); 

forT =1:ntra-1 


Qi ene: 


res 


ui(j) = ui(j) +0.5 * (DF(T +1) * V,4(T +1) + Df (T) * V,i(T)) * abs(V,i(T +1) — V,i(T))/(ni(y) +. 


end 


clearV,iT 
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end 


(g) Main command file for DFA 


close all, clear all, global niL vmLi vtfi kB Tfi Mi vmaxLi vcLi Ef dx X nx alpha 


% input(‘Type the name of your input file:’); 

Input; 

CouplingPSandS; % gives input parameters for sheath simulation, 
uLe = —vt fe «x De/Ce/sqrt(pi); uLi = umLi — vt fi * Di/Ci/sqrt(pi); 
DLe = sqrt(epsilon « kB « T fe/neL/e?); 

L =eval(L); 

vmazr Li = Taumin * vt fi — umL1; 

I =1,; % Initializing the number of iteration counter 

X = linspace(0, L,nx); 

dx = X(2) — X(1); 

V Li = linspace(umazxLi, vcLi, ntra); 

dV Li = abs(V Li(2) — VLi(1)); 

ifdz/DLe > 8 

disp(‘ Attention: Your grid width is comparable/larger than Debye length.’) 
dx , DLe 

disp(‘Press ENTER to continue or CTRL-C to stop.’) 


pause 
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end 

phi = eval(phi0z); phi(1) = phi0; phi(nx) = phiL; % load potential 
forj=2:nz-1 

Ef(j) = .5* (phi(j — 1) — phiG + 1))/da; 

end 

Ef (1) = 5 « (3 * phi(1) — 4 * phi(2) + phi(3))/dz; 

Ef (nz) = .5 * (—3 * phi(nx) + 4 * phi(nx — 1) — phi(na — 2))/dz; 

DF Aiondensity 

ne = neL,exp(e * phi/kB/T fe). * (1+ erf(sqrt(e * (phi — phi(1))/kB/T fe))); 


ni = niL « ni/max(ni); ne = niL * ne/maz(ne); nin = ni/niL; ui = uiLl * ui/max(ut); rho = Ze * ni - 


d = —2 x ones(1, nz); d(1) = 1;d(nx) = 1; u = ones(1, nx — 1); 

ML = diag(d) + diag(u, 1) + diag(u, —1); ML(1, 2) = 0; ML(nz, nx — 1) = 0; 
clear d u 

Poissonsolver; while max(abs(phin-phi)) > deltaphi 

fluctuation = max(abs(phin — phi)) 

phi = phin; 

fore =2 ne 1 

Ef(j) = .5* (phi(j — 1) — phiG + 1))/da; 


end 
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Ef (1) = 5 * (3 * phi(1) — 4 * phi(2) + phi(3))/da; 

Ef (nx) = .5 * (—3 * phi(nx) + 4 * phi(nx — 1) — phi(nax — 2))/da; 

DF Aiondensity; 

ne = neL,exp(e * phi/kB/T fe). * (1+ erf(sqrt(e * (phi — phi(1))/kB/T fe))); 
ni = niL * ni/maz(ni); ne = niL * ne/maz(ne); 

nin = ni/max(nt); 

nen = ne/max(ne); 

ucLe = sqrt(—2 « e * phi0/Me); 

ui = uli * ui/mazx(ut); 

rho= Zexni—ex ne; 

poissonsolver; 

end 

phi = phin; 

DF Aiondensity; 

ne = neL,exp(e * phi/kB/T fe). * (1+ erf(sqrt(e * (phi — phi(1))/kB/T fe))); 
ni = nil *« ni/max(ni);ne = nil x ne/max(ne); nin = ni/niL; nen = ne/max(ne); 
ui = uiL *« ui/max(ut); uin = ui/max(ut); 

rho = Zexni—exne; 

x = X/DLe: 

Ef =Ef,Efn=Ef/max(Ef); 


plot(x, phin, ‘r’), xlabel(‘normalised distance’), ylabel(‘phin’) 
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